!p.charsize=1.5

l=188
rmin=0.6
rmax=0.96
re=rmin+(rmax-rmin)*indgen(l)/double(l-1)
Rg=13732.
cp=2.5*rg
sr=6.95989e8

filename='sms.dat'
dd=read_ascii('sms.dat')
r=dd.field1[0,*]
c=dd.field1[1,*]
rho=dd.field1[2,*]
p=dd.field1[3,*]
gamma=dd.field1[4,*]
t=dd.field1[5,*]

;only between 0.6 and 0.97

ii=where(r ge rmin and r le rmax)
r=reverse(r[ii])
c=reverse(c[ii])
rho=reverse(rho[ii])
p=reverse(p[ii])
gamma=reverse(gamma[ii])
t=reverse(t[ii])

c_new=interpol(c,r,re)
rho_new=1000.*interpol(rho,r,re)
p_new=interpol(p,r,re)/10.
gamma_new=interpol(gamma,r,re)
t_new=interpol(t,r,re)

print,rho_new[0]
print,t_new[0]
print,p_new[0]


s_new=alog(p_new)/gamma_new-alog(rho_new)
th=exp(s_new/cp)

jj=where(re ge 0.7 and re le 0.71)
th00=mean(th[jj])
t00=mean(t_new[jj])

;th_new=t_new*(p_new[0]/p_new)^((gamma_new-1.)/gamma_new)
th_new=t_new*(p_new[0]/p_new)^(2./5.)

plot,re,th_new,yr=[min(th_new),max(th_new)],xr=[0.6,1]
oplot,[rmin,rmax],[t00,t00],li=2
openw,1,'strat.dat'
for i=0,l-1 do printf,1,rho_new[i],t_new[i],p_new[i],th_new[i]
close,1
save,filename='solardata.sav',re,rho_new,t_new,p_new

end
